library(tidyverse)
library(readstata13)
library(foreign)
library(readxl)
library(writexl)
library(lubridate)
library(PanelMatch)

rm(list=setdiff(ls(), "data_final"))

setwd("/Volumes/Backup Plus/Storage-OrganizedFiles/Article1-ReplicationMaterials_FullReplication/3.RobustnessChecks")
getwd()

data1 <- read.dta13("Ortega(2024)-DataAnalysis_Updated(09302024).dta")

data_final <- data_final %>%
  filter(year > 1987)

#################################################################################
# Analysis
#################################################################################

# IV: Protests
# A: Create the Matched Set (M)

pm.none5.p <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv2.p_strat", refinement.method = "none",
                         data = data_final, match.missing = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.2.no_p_str, 1:3)) 
                         + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv2.1.kciv_FARC_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE,
                         use.diagonal.variance.matrix = TRUE)

pm.maha5b.p <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv2.p_strat", refinement.method = "mahalanobis",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.2.no_p_str, 1:3)) 
                          + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                          size.match = 5, qoi = "att" , outcome.var = "dv2.1.kciv_FARC_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

pm.ps.m5.p <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv2.p_strat", refinement.method = "ps.match",
                         data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                         + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.2.no_p_str, 1:3)) 
                         + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv2.1.kciv_FARC_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE)

pm.ps.w5.p <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv2.p_strat", refinement.method = "ps.weight",
                         data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                         + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.2.no_p_str, 1:3)) 
                         + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv2.1.kciv_FARC_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE)

pm.CBPSm5.p <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv2.p_strat", refinement.method = "CBPS.match",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.2.no_p_str, 1:3)) 
                          + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                          size.match = 5, qoi = "att" , outcome.var = "dv2.1.kciv_FARC_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

pm.CBPSw5.p <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv2.p_strat", refinement.method = "CBPS.weight",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                          + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.2.no_p_str, 1:3)) 
                          + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                          size.match = 5, qoi = "att", outcome.var = "dv2.1.kciv_FARC_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE)

# IV: Violent Self-Protection
# A: Create the Matched Set (M)
pm.none5.v <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv4b.v_strat", refinement.method = "none",
                         data = data_final, match.missing = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.4b.no_v_str, 1:3)) 
                         + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv2.1.kciv_FARC_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE,
                         use.diagonal.variance.matrix = TRUE)

pm.maha5b.v <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv4b.v_strat", refinement.method = "mahalanobis",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.4b.no_v_str, 1:3)) 
                          + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                          size.match = 5, qoi = "att" , outcome.var = "dv2.1.kciv_FARC_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

pm.ps.m5.v <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv4b.v_strat", refinement.method = "ps.match",
                         data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                         + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.4b.no_v_str, 1:3)) 
                         + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv2.1.kciv_FARC_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE)

pm.ps.w5.v <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv4b.v_strat", refinement.method = "ps.weight",
                         data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                         + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.4b.no_v_str, 1:3)) 
                         + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv2.1.kciv_FARC_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE)

pm.CBPSm5.v <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv4b.v_strat", refinement.method = "CBPS.match",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.4b.no_v_str, 1:3)) 
                          + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                          size.match = 5, qoi = "att" , outcome.var = "dv2.1.kciv_FARC_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

pm.CBPSw5.v <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv4b.v_strat", refinement.method = "CBPS.weight",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                          + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.4b.no_v_str, 1:3)) 
                          + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                          size.match = 5, qoi = "att", outcome.var = "dv2.1.kciv_FARC_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE)

# IV: Sanctuary
# A: Create the Matched Set (M)
pm.none5.s <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv3.s_strat", refinement.method = "none",
                         data = data_final, match.missing = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.3.no_s_str, 1:3)) 
                         + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv2.1.kciv_FARC_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE,
                         use.diagonal.variance.matrix = TRUE)

pm.maha5b.s <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv3.s_strat", refinement.method = "mahalanobis",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.3.no_s_str, 1:3)) 
                          + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                          size.match = 5, qoi = "att" , outcome.var = "dv2.1.kciv_FARC_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

pm.ps.m5.s <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv3.s_strat", refinement.method = "ps.match",
                         data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                         + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.3.no_s_str, 1:3)) 
                         + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv2.1.kciv_FARC_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE)

pm.ps.w5.s <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv3.s_strat", refinement.method = "ps.weight",
                         data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                         + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.3.no_s_str, 1:3)) 
                         + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv2.1.kciv_FARC_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE)

pm.CBPSm5.s <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv3.s_strat", refinement.method = "CBPS.match",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.3.no_s_str, 1:3)) 
                          + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                          size.match = 5, qoi = "att" , outcome.var = "dv2.1.kciv_FARC_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

pm.CBPSw5.s <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv3.s_strat", refinement.method = "CBPS.weight",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                          + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.3.no_s_str, 1:3)) 
                          + I(lag(dv2.1.kciv_FARC_r, 1:3)),
                          size.match = 5, qoi = "att", outcome.var = "dv2.1.kciv_FARC_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE)

#### Plots ####

msets.none5.p <- pm.none5.p$att
msets.maha5b.p <- pm.maha5b.p$att
msets.ps.m5.p <- pm.ps.m5.p$att
msets.ps.w5.p <- pm.ps.w5.p$att
msets.CBPSm5.p <- pm.CBPSm5.p$att
msets.CBPSw5.p <- pm.CBPSw5.p$att

# C Panel Estimate
pe.results5.p2_maha5b <- PanelEstimate(sets = pm.maha5b.p, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

pe.results5.p2_ps.m5 <- PanelEstimate(sets = pm.ps.m5.p, data = data_final,
                                      se.method = "conditional",
                                      number.iterations = 1000,
                                      confidence.level = .95)

pe.results5.p2_ps.w5 <- PanelEstimate(sets = pm.ps.w5.p, data = data_final,
                                      se.method = "conditional",
                                      number.iterations = 1000,
                                      confidence.level = .95)

pe.results5.p2_CBPSm5 <- PanelEstimate(sets = pm.CBPSm5.p, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

pe.results5.p2_CBPSw5 <- PanelEstimate(sets = pm.CBPSw5.p, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

msets.none5.v <- pm.none5.v$att
msets.maha5b.v <- pm.maha5b.v$att
msets.ps.m5.v <- pm.ps.m5.v$att
msets.ps.w5.v <- pm.ps.w5.v$att
msets.CBPSm5.v <- pm.CBPSm5.v$att
msets.CBPSw5.v <- pm.CBPSw5.v$att

pe.results5.v2_maha5b <- PanelEstimate(sets = pm.maha5b.v, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

pe.results5.v2_ps.m5 <- PanelEstimate(sets = pm.ps.m5.v, data = data_final,
                                      se.method = "conditional",
                                      number.iterations = 1000,
                                      confidence.level = .95)

pe.results5.v2_ps.w5 <- PanelEstimate(sets = pm.ps.w5.v, data = data_final,
                                      se.method = "conditional",
                                      number.iterations = 1000,
                                      confidence.level = .95)

pe.results5.v2_CBPSm5 <- PanelEstimate(sets = pm.CBPSm5.v, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

pe.results5.v2_CBPSw5 <- PanelEstimate(sets = pm.CBPSw5.v, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

msets.none5.s <- pm.none5.s$att
msets.maha5b.s <- pm.maha5b.s$att
msets.ps.m5.s <- pm.ps.m5.s$att
msets.ps.w5.s <- pm.ps.w5.s$att
msets.CBPSm5.s <- pm.CBPSm5.s$att
msets.CBPSw5.s <- pm.CBPSw5.s$att

pe.results5.s2_maha5b <- PanelEstimate(sets = pm.maha5b.s, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

pe.results5.s2_ps.m5 <- PanelEstimate(sets = pm.ps.m5.s, data = data_final,
                                      se.method = "conditional",
                                      number.iterations = 1000,
                                      confidence.level = .95)

pe.results5.s2_ps.w5 <- PanelEstimate(sets = pm.ps.w5.s, data = data_final,
                                      se.method = "conditional",
                                      number.iterations = 1000,
                                      confidence.level = .95)

pe.results5.s2_CBPSm5 <- PanelEstimate(sets = pm.CBPSm5.s, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

pe.results5.s2_CBPSw5 <- PanelEstimate(sets = pm.CBPSw5.s, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

save.image(file='RCh3-ODDHP_AlternativeModel_FARC.RData')

########################################################
## End of this part
########################################################